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Background, foreground and nearby matter influence on strong 
gravitational lenses 
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ABSTRACT 

We investigate strong lensing by non-singular finite isothermal ellipsoids taking into account 
the influence of the matter along the line of sight and in the close lens vicinity. We com- 
pare three descriptions of light propagation: the full approach taking into account all matter 
inhomogeneities along the rays, the single plane approach, where we take into account the in- 
fluence of the strong lens neighbours but neglect the foreground and background objects, and 
the single lens approach. In each case we simulate many strong lensing configurations placing 
a point source at the same redshift but in different locations inside the region surrounded by 
caustics. We further analyze configurations of four or five images. For every simulated strong 
lensing configuration we attempt to fit a simplified lens model using a single isothermal ellip- 
soid or a single isothermal ellipsoid with external shear. The single lens fits to configurations 
obtained in the full approach are rejected in majority of cases with 95% significance. For 
configurations obtained in the single plane approach the rejection rate is substantially lower. 
Also the inclusion of external shear in simplified modeling improves the chances of obtaining 
acceptable fits, but the problem is not solved completely. The quantitative estimates of the 
rates of rejection of simplified models depend on the required accuracy of the models, and we 
present few illustrative examples, which show that both matter close to the lens and matter 
along the rays do have important influence on lens modeling. We also estimate the typical 
value of the external shear and compare the fitted parameters of the simplified models with 
the parameters of the lenses used in the simulations. 

Key words: gravitational lensing: strong and weak - large-scale structure of the Universe 



1 INTRODUCTION 

Strong gravitational lensing of distant sources has low probability 
and is usually attributed to a single concentration of matter close 
to the line of si ght, while the in fluence of anything else is ne- 
glected (see e.g. Kochanek 200GS for a review of strong lensing). 
In some cases simple, single component models of lenses are not 
satis factory. Even the fi r st lensed quasar o bserved, QSO 0957+561 
A,B ( Walsh et alj[l979t iGunn et al.lll979h requires a combination 
of a galaxy and a galaxy cluster to model its basic properties 
dYoungetalJll98ih . 

The influence of another object along the line of 
sight on the main lens prop er ties has bee n inve stigated by 
iKochanek & Apostolaki3 dl988l) . iKeeton et all tl99l) argue that 
adding an external shear to an elliptic lens model greatly improves 
the fits of four image configurations, but th ey conclude tha t 
the shear is probably related to the main lens. iBar-Kanal dl99fj|) 
investigates the influence of the density fluctuations along the 
line of sight on the accuracy of Hubble constant determination 
based on time delays measured in some of the multi-image QSO 
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systems. Ichen et all d2003h check the influence of the substructure 
not rela ted to the main len s on th e observed flux ratios in lensed 
images. Wambs ganss et all d2004) study the probability of giant 
arcs using the results of a large cosmological simulation and 
multi-layer approach to backward ray tracing. They find that 
the matter on the way sli ghtly increases the probability of arcs. 
IWambsganss et all d2005h address specifically the question of 
secondary lenses on the line of sight. According to this study the 
role of secondary lenses is a strong function of source redshift: 
it can be neglected in 95% of cases for a source at zs = 1, but is 
important in 38% of cases at zs = 3.8. Using single plane approach 
introduces system atic errors to the estimated l ens parameters. 
In a recent paper iD'Aloisio & Nataraianl d201ll) investigate the 
influence of the secondary lenses in the line of sight on the 
accuracy of estimating the dark energy density. 

It is probably impossible to quantify the influence of the mat- 
ter along the line of sight and of the strong lens close (but possibly 
unobservable) neighbours in a single case of real observed multiple 
image configuration. In this paper we are going to obtain a statisti- 
cal measure of such an influence performing simulations of many 
multiple image configurations using a realistic model of light prop- 
agation in an inhomogeneous Universe model. To find the role of 
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matter not belonging to the strong lens, we follow the rays includ- 
ing or not the lens neighbours and inhomogeneities along the rays. 
Again, the direct comparison between individual image configura- 
tions obtained as a result of different descriptions of light propa- 
gation seems useless, since one should compare positions and flux 
ratios of images of a point source at the same position and it is not 
clear how to compare source positions between the different mod- 
els of light propagation. If image configuration is considered, it is 
the source position relative to the caustic structure which counts 
and since the caustic structures are different the problem persists. 

Instead of making comparisons on a one to one basis, we com- 
pare the properties of all four or five image configurations between 
the models of light propagation. If the influence of matter outside 
the strong lens were unimportant, each lens - image configuration 
would be similar to some configuration obtained with a single lens 
model. To check whether this is the case we attempt to fit such 
simplified model to all image configurations considered. The suc- 
cess rate of the fitting procedure applied to configurations obtained 
with a different model of light propagation statistically measures 
the similarity of this model to the single lens model. 

In Sec. 2 we describe our models of light propagation. Sec. 3 
presents tools used to compare different models and the results of 
such comparison. Discussion and conclusions follow in Sec. 4. 



2 MODEL OF THE LIGHT PROPAGATION 
2.1 Deflection angles 

The results of the Millennium Simulation dSpringel et al.ll2005h 
give the matter distribution (the positions and velocities of ~ 10 10 
simulation particles) in the 64 epochs corresponding to the redshifts 
Zi ^ 127. In our calculations we assume, that the propagation of 
a ray between redshifts z,-, = (zj_i + z/)/2 and Zi 2 = (z;+i +zi)/2 is 
affected by matter distribution given for the epoch z,- We are inter- 
ested in sources at high redshifts (z Js 1), so there are always many 
such matter layers on the ray path, and e ach of them can be consid- 
ered thin. Thus the multilayer approach (Schneider & Weiss ll988l) 
with layers corresponding to the Millennium epochs is natural. The 
position of the ray in the N-th layer in the angular coordinates is 
given as 



Pw = Pi 



N-l 

I 

/=1 



«f(Pi) 



(1) 



where djN is the angular diameter distance as measured by an ob- 
server at epoch i to the source at epoch N, dfj - the angular diameter 
distance to the same source measured by a present (z = 0) observer, 
and 0t/(p,) is the deflection angle in the i-th layer at the position 
P ; . In a flat cosmological model the angular diameter distances in 
the lens equation can be replaced by comoving distances. Below 
we use comoving distances D(z), denoting £>, = D(zi) for short. 
Since the comoving distances are additive in a flat model, one has 
diN/dN = (Dn —Di)/Dn- I n the calculations we ap ply more effi- 
cient recurrent formula of ISeitz & Schneider! dl992b . equivalent to 
the above equation. 

To calculate the deflection angle in a given layer we need a 
description of its matter distribution. This is done in two steps. 
The averaged matter density is defined on a coarse grid of 256 3 
cells. The gravitatio nally bound haloes have been described by 
|Pe Lucia & Blaizotl d2007l) and iBertone et all d2007h . Both kinds 
of data are accessible from the Virgo - Millennium Da tabase 
dLemson &" Springel 2006). Ijaroszvns ki & Kostrzewai d2010h use 



the same information on matter distribution, but their calculation 
of lensing effect is based on different approach. 

Since only the positions, virial masses, virial radii, and virial 
velocities of the haloes can be obtained from the Database, the de- 
tailed matter distributions of each halo must be supplemented. We 
use non singular isothermal ellipsoids (NIE) as models of individ- 
ual lenses. The averaged matter density includes also the matter 
belonging to haloes, so the halo mass models must be compensated 

by a shallow negative density distribution (see below); 

In constructing the matter layers we follow ICarbone et al.l 

(2008), randomly rotating and shifting simulation cubes corre- 
sponding to different epochs, which eliminates the consequences 
of periodic boundary conditions used in their calculations. The 
comoving thickness of the i-th layer £),, — Dj t is of the order of 
10 2 Mpc, always smaller than the simulatio n cube size 500//i M pc 



The 2D layers are not periodic in general. iHilbert et al.l d2009l) de- 
fine specific projection directions to obtain periodic 2D density dis- 
tributions and use 2D Fourier transforms in calculations of the de- 
flection angles. Our 3D grid (256 3 ) is small enough to employ spec- 
tral methods. Using Poisson equation we calculate 3D gravitational 
acceleration g on the grid. The deflection of a ray passing the i- 
th layer at position p, due to the averaged (or background) matter 
distribution is given as: 



"■ CS (P,0 



.ID:, 



g± (D,fr)dD 



(2) 



where g ± is the component of gravitational acceleration perpen- 
dicular to the ray, and its value at any location is obtained by the 
interpolation on the 3D grid. (D, P) are used as coordinates. Since 
our simulated maps of the sky cover typically regions of a few min- 
utes of arc in size, the position components (D$ x ,D$y) may be 
treated as Cartesian coordinates in the plane perpendicular to the 
propagation direction, and the comoving distance D serves as ra- 
dial coordinate. Because the integration should be performed over 
the proper distance, instead of comoving one, the result is divided 
by the factor 1 + Zi ■ 

For each individual halo we a pply the model of non -singular 
isothermal ellipsoid as described bv lKormann et al. dl994l) . The de- 
flections in 2D real notation are given bv lKochaneM d2006l) : 
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where: 



a>(x,y,q,r ) 



q 2 (x 2 + rZ)+y 2 q' 



q- 



The ray crosses the lens plane at (x.y), the lens centre is placed 
at the origin of the coordinate system, the major axis along x. The 
axis ratio is given by q, rn is the core radius, and ocq is the deflection 
angle parameter (see below). 

The convergence K (which is proportional to the surface mass 
density) has the required (NIE) form: 

1 



K{x,y,q,r ) = -(a^ + Oyj) 



Ok) 



2 u>(x,y,q,r ) 



(4) 



Far from the centre the 2D matter distribution becomes elliptical, 
as required. 

The NIE matter distribution is infini te with divergen t total 
mass. To obtain a finite mass lens we follow lKochaneS 1 2006) com- 
bining two NIE distributions with different core radii r\ Cr2. In- 
side r\ the surface density is approximately constant, isothermal in 
between, and sharply falling outside r^'- 
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Figure 1. The deflection angle of the compensated lens. For illustration 
we use: r\ = 1, r% = 10 2 n, and /-/;„, = 10 2 r2. The plots show the radial 
deflection component ct,- as function of radius for fixed polar angle <S) = 90°, 
60°, 30°, and 0° (top - down). Dashed lines are used for q = 0.3, and solid 
lines for q — 0.7. 



K 0, y, i ) = K 0, y, q, n ) - k(x, y, q, r 2 ) 



1 (Xo(^-^) 

4 O) 3 (x,j,^,0) 



(5) 



which gives a non-singular, finite mass lens with the deflection an- 
gle: 



<*(x,y,q) =a(x,y,q,n)-a(x,y,q,r 2 )->-a (r2-ri)- 



(6) 



For large radii r 2> f2 the lens acts as a spherical finite mass con- 
tained within r 2 . We identify the mass and radius of the lens with 
the virial mass and radius of the halo, so the characteristic deflec- 
tion angle is given as: 



OCq 



4GM V 



4GM V 
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(7) 
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To compensate the lens we use a circular disk of radius rr im 3> 
r v ; r with constant surface density and with the deflection angle: 
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4GM V 
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(8) 



where 

4 ^ 

defines the size of the disk, pg is the average matter density in the 
Universe at the lens redshift. In our approach we always use r\ = 
0.01r2, so we avoid singularity at the lens centre, but we do not 
assign any physical meaning to this parameter. 

Finally we check numerically the dependence of the deflection 
angle on various parameters. We use the formula for a compensated 
lens: 



O(r) = a{x,y,q,n)-a(x,y,q,r 2 )-ac m Jr 



(10) 



We illustrate the dependence of deflection angle on position in 
Fig. U] using arbitrary units and a case with r\ <C r 2 <C Him- The 



compensation is not exact, but at r ~ r\[ m the deflection angle drops 
by two orders of magnitude. 

A compensated lens has only a finite range. Its radius depends 
on the properties of each individual halo, but investigating the de- 
flection by haloes at any position in a given layer we take into ac- 
count only the finite number of lenses: 



a haloes (p) = £ o(p p , } where . D .|p_ p ^ r li 



Jim 



(ii) 



The expression 0C(($ — fij) stands for the deflection at the position 
P caused by the lens at position p ; - within its limiting range. 
The total deflection in the i-th layer 



a ! (p i ) = af g (P i ) + al ialoes (P ; ) 

can be calculated for a ray passing at any location p ; . 



(12) 



2.2 Maps of the sky 

We calculate deflection angles for all layers of interest storing them 
on the2048 2 grids with 0.1 arcsec resolution. Using a set of stored 
deflection angles for all layers and interpolating on the grids one 
can apply Eq.(l), mapping a region of ~ 3' x 3' in the sky. We 
investigate light propagation from sources at zs ~ 2 taking into ac- 
count the influence of several thousands haloes in the field of view 
and some in adjacent regions, if they are within their ri; m ranges. 
We construct eight separate maps pointing in random directions. 

When investigating multiple image properties we need 
zoomed maps of smaller parts of the sky. For better resolution we 
use finer grids with deflections interpolated from previous calcula- 
tions with the help of bi-cubic spline, so the interpolated deflection 
derivatives are continuous. The backward ray shooting on a finer 
grid is repeated. 

Treating each of the individual haloes in the field of view as 
an isolated lens of known velocity dispersion and using SIS model 
we calculate its Einstein radius r£ for a source at given redshift zs 
in a homogeneous Universe model. In each simulated map of the 
sky we find 10 lenses with the largest Einstein radii. We expect that 
these dominating lenses have the highest probability of producing 
multiple images. 

In the single SIS lens case the second image is possible only 
for a source within the Einstein ring. The brighter of the two images 
lies within 2 r£ from the lens centre and the dimmer within 1 rg, 
so we expect that in more complicated problem, including lens el- 
lipticity and the influence of other lenses the multiple images lie 
within few Einstein radii from each other. 

We check for the presence of other haloes inside the circle of 
the radius 3 ye surrounding a dominating lens. If they are present 
we enlarge the region of interest including 3 ye zones around all 
companions. Finally we repeat backward ray shooting inside a 
square on the sky including the region of interest. The fine grids 
giving the deflection angles in consecutive layers encompass still 
larger areas, so one can follow rays deflected off the main region. 
The result of the ray shooting is a vector array: 



Vn 



MP! 



(13) 



where f$ gives the positions in the source plane of rays apparently 
coming from the directions Pj on the observer's sky. Superscripts 
k, I enumerate the rays. The Pn(Pi ) relation is given in Eq.(l). 
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Figure 2. Two examples of critical lines and caustic structures obtained using Case 1 light propagation model. Critical lines are shown in the left column. 
The coordinates are relative to the centre of 3' x 3' maps and are expressed in arcseconds. The corresponding caustic structures are shown on the right. The 
perturbations of the caustic structure by the surrounding haloes are apparent in the example shown in the upper row. In the lower row the action of two close 
elliptical lenses of similar masses, having different position angles completely changes the shapes of critical lines and caustics. 



2.3 Three propagation models 

To estimate the influence of the foreground and background matter 
in the beam and also of the lens neighbours belonging to the same 
redshift layer we consider three scenarios f or the backward ray 
shooti ng. Very similar method is employed bv lPuchwein & Hilbertl 
(2009). In the full approach (Case 1) we include deflection of rays 
by matter in all layers between the source and observer, as de- 
scribed above. In the one layer approach (Case 2) we include de- 
flection by all matter in the layer of dominating lens only, neglect- 



ing the influence of other layers. In the single lens approach (Case 
3) only the deflection by the dominant lens is taken into account. 



2.4 Image finding 

To find critical lines in the sky we solve the standard equation 
JSchneideretal.ll992l) : 



det A = 



A = 



9Pa 



(14) 
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where A is the deformation matrix. The equation is solved approxi- 
mately on the grid and then refined by iterations. Finding of caustic 
lines in the source plane is then straightforward. 

Usually the critical lines and caustic structures obtained in 
three Cases of backward ray shooting are similar and the influence 
of matter outside the dominating lens is not apparent. There are 
however several cases showing the influence of smaller mass haloes 
on the critical lines and (to a lesser extent) caustics as shown in the 
upper row of Fig. [2] An example of a strong caustic perturbation is 
shown in the lower row of Fig. [2] 

We start image finding using approximate methods on the grid. 
We replace the point source by an extended surface brightness pro- 
file with the Gaussian shape and the characteristic radius of few 
pixels. The related surface luminosity in the sky is given by: 

/ obs (P 1 )=/ src (P w (P 1 )) (15) 

The local maxima of the observed surface luminosity are the posi- 
tions of the images of the source centre. We find brightness maxima 
on the grid and use them as approximate solutions to the lens equa- 
tion. The improved positions are obtained by iterations. 

Using the above method we find multiple images of point 
sources placed behind lenses. For each lens we try few hundreds 
to several thousands source positions (depending on the size of 
the caustics region) distributing them evenly. Each position consid- 
ered represents the same surface area in the source plane. Finding 
a source at any location has the same probability if they are dis- 
tributed randomly. (We neglect here the amplification bias caused 
by the strong lensing.) The procedure of image finding for many 
source positions is repeated for each of the ten strongest lenses be- 
longing to each of the eight simulated maps of the small fragments 
of the sky obtained with the three different models of light propa- 
gation. 



3 FITTING THE SIMULATED STRONG LENSES WITH 
SIMPLIFIED MODELS 

In the case of pure elliptical, non-singular lens one would expect 
(in general) 1, 3, or 5 images. Our numerical method using at the 
beginning a small extended source and having finite resolution may 
miss some of the images. (See the Discussion). In Cases 1,2 the 
external influence may also change the number of images. Since 
configurations with more images probably better constrain the lens, 
we concentrate on cases with five or four images. 

To estimate the influence of the galaxies in the lens vicinity 
and the matter along the line of sight on the properties of the strong 
lensing, we attempt to model all our simulated cases of multiple 
imaging using a single lens model in a uniform Universe. The sin- 
gle lens we are using in modeling is a non-singular finite isothermal 
ellipsoid used also in the simulations for each of the haloes. When 
tracing the rays in simulations we interpolate all deflection angles 
from earlier computed arrays. Single lens modeling uses analyti- 
cally calculated deflection angles and their derivatives, so it serves 
also as a test of ray tracing simulations in Case 3. 

We also try a more sophisticated lens model using the same 
non-singular finite isothermal ellipsoids with external shear. In this 
way the tidal influence of masses close to the rays is at least par- 
tially represented. 

Simulations give the lens position, image positions and image 
amplifications corresponding to their energy fluxes. These param- 
eters would also be observed in reality. The source position P s is 
unknown but it must be the same for all the images, and we use it as 



a model parameter. The source luminosity is also unknown, but the 
observed fluxes must be proportional to the lens amplifications in 
corresponding positions. Thus the model has to reproduce ampli- 
fication ratios and not their absolute values. The intrinsic lens pa- 
rameters (axis ratio q, characteristic deflection angle oio, the virial 
radius r v j r , and the lens position angle in the sky <])) are not so easy 
to measure and we treat them as unknown model parameters. Also 
the shear components yi , j2 belong to the class of intrinsic param- 
eters. Below we use the value of shear defined as 



(16) 



We treat the lens position, image positions and image fluxes 
ratios given by the simulations as observed quantities, each with 
its own uncertainty. We attempt to fit every simulated strong lens- 
ing case with our simplified model. We are going to reproduce the 
observed lens position, image positions and image flux ratios (or 
differences in image stellar magnitudes) looking for a minimum of 
X 2 taking all three constraints into account: 



X 2 — Xl + Xi + Xm 



(17) 



The first term, controlling the lens position, used in the model has 
the obvious form: 

x 2 = (Pl-Pl°) 2 (18) 

where the subscript L stands for "lens" and Ol is the assumed ac- 
curacy of the lens position. 

Using the simplified model we can calculate the source posi- 
tions p s ' related to observed image positions P T '°. The deformation 
matrix Am and magnification matrix M/s can be calculated at each 
image position: 



Ps = Ps(Pi 



A (0 = 



(19) 



The mismatch between pg and P s implies magnified mismatch be- 
tween modeled a nd observed image positions, which gives for %} 
dKochaneldbood) : 

I M (i r (P^-Ps]' 2 



The fitting statistic for flux ratios is given as: 
. . im' — m"' — <m' — ///' 

Xm = 2j ZJ2 



(20) 



(21) 



In all cases the quantities with extra superscript "0" are taken from 
simulations and mimic the observed values while the quantities 
without this superscript are given by the model. 

We use the following characteristic values for the accuracy 
parameters: = 0.01 arcsec, Ol = 0.1 arcsec, and o" m = 0.1 or 10 
(in stellar magnitudes) This choice takes into account the fact that 
image positions are easier to measure than the lens position. With 
the small value of the flux measurement error (a m = 0.1) we use 
the fluxes as important factor in modeling. The other value makes 
the fluxes unimportant in the model fitting. Using it one can check 
to which extent models well reproducing image positions are also 
good in reproducing ratios of their fluxes. 

We start fitting process using a model with all parameter val- 
ues, also the intrinsic lens parameters, taken from the simulation. 
This is not a general approach, but it is fast in finding solutions 
close to the starting point in parameter space (if they are present). 
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Table 1. Acceptability of fits - dependence on the model 



Table 2. Dependence of the acceptability of fits on required model accuracy 



Shear no no yes yes 

a m 0.1 10. 0.1 10. 

Casel 0.08±0.13 0.11 ±0.13 0.43±0.17 0.68±0.22 

Case 2 0.38±0.21 0.41 ±0.18 0.59±0.22 0.78±0.19 

Case 3 .9998 ±.0002 .9998 ±.0001 .9997 ±.0002 .9998 ±.0001 



Note: The table shows the dependence of the rate of acceptability of fits 
on the method of treating the light propagation (Cases 1-3, see the text 
for details), for models neglecting ("no") or taking into account ("yes") the 
external shear and reproducing (o~ m = 0.1) or not (o~ m = 10) the flux ratios. 
The assumed errors in fitting the positions of the lens center and of the 
images are kept constant (o~l = 0.1, o~i = 0.01). 

It is also our goal to obtain the "observable" parameters as close to 
simulated values as possible, but the intrinsic lens parameters are 
not constrained a' priori and have no preferred values. 

Our approach to model fitting is not general, since we are 
probing only the small region of parameter space and there is a 
chance that we miss some acceptable models of the lens. On the 
other hand, if the influence of the environment were weak, the re- 
sulting lens model should be close to the original lens used in the 
simulations. Failure to find a model which is similar to the original 
can be inteipreted as the result of the strong influence of the matter 
in the rays vicinity. 

The number of the degrees of freedom (DOF) depends on the 
number of images (N lm ) and on the number of intrinsic model pa- 
rameters, which is iV par = 4 for an elliptic lens (q, Oq, r wu , <])) and 
Nnai = 6 if the external shear (ji, Y2) is taken into account. The 
lens position is both observed and modeled, so it has no impact. The 
source position requires 2 parameters, and the source internal lumi- 
nosity (which we do not use explicitly in modeling) - another one. 
When modeling the flux ratios with a m = 0.1 we effectively use the 
magnification information from the simulations, but for o m = 10 it 
is practically lost. In calculations we use: DOF = 3(A , ; m — 1) — iVpar 
(if <T m = 0.1) or DOF = 2{N im - 1) -2V par (if <T m = 10). (That im- 
plies that modeling of four-image configurations including shear 
and neglecting flux ratios is impossible, since DOF — in this 
case). Using tables of % 2 distribution we reject fits if they give 

X 2 >xl^{DOF) (22) 

where %q 95 corresponds to 95% significance for a given DOF. 

The influence of the matter distribution outside the lens on 
the lensing process can be characterized by the difficulty in ob- 
taining an acceptable single lens model in such more complicated 
situations. To measure the rate of success in fitting procedure we 
assign the same weight to each investigated source position since 
they are distributed evenly in the source plane. Thus the more mas- 
sive lenses have greater influence on the results as compared with 
less massive. This is the simplest consequence of the observational 
selection. More subtle effects, like magnification bias, influence of 
image configuration etc are too difficult to model within our ap- 
proach. 

Since we perform our simulations in eight different, randomly 
chosen regions of the synthetic Universe, we can estimate the role 
of the cosmic variance in our simulations. Calculating the rate of 
success for each of the regions separately and finding the dispersion 
of the results, we get an error estimate. 

In Table. 1 we show the chances that a simulated strong lens- 
ing case has an acceptable simplified model. In Case 3 acceptable 
models are possible to obtain in almost 100% of simulated lens- 
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Note: The table shows the dependence of the rate of acceptability of fits 
on the method of treating the light propagation (Cases 1-3, see the text 
for details) and on the assumed errors in measured positions of the center 
of the lens (o~l) and of the images (o~i). All models reproduce flux ratios 
(rj m = 0.1) 



image configurations, which shows that our fully numerical (used 
in obtaining lens - image configurations) and semi-analytical (used 
in fitting models) descriptions of the same physical approach are in 
complete agreement. 

The results depend on the adopted values of position accu- 
racy, (see below) but the difficulty in modeling realistic configura- 
tions using a single elliptical lens without external shear is appar- 
ent, even if the reproduction of correct flux ratios is not required. 
The inclusion of external shear substantially improves the chances 
of obtaining acceptable solutions with simplified lens models, but 
still there is no guarantee of success, especially if one requires the 
reproduction of the observed flux ratios. 

The present day astrometry gives better accuracy of image po- 
sitions than assumed in majority of our simulations. On the other 
hand the micro-lensing may change the measured fluxes of individ- 
ual images in a random way, which diminishes their role as model 
constraints. In Table. 2 we show the results of our approach with 
increased required accuracy of positions of the lens, the images, or 
both, for models with shear, reproducing flux ratios. 

In the following Fig. [3] we show the probability P(> y) that 
an acceptable fit uses external shear value greater than Y- The plots 
refer to fits taking into account image positions but not their flux 
ratios. 

The influence of lens neighbours and matter in the lens fore- 
ground/background is also illustrated by the changes in the intrinsic 
lens parameters when comparing the fitted model with the original 
lens used in the simulation. In reality one has only one set of lens 
parameters obtained from modeling and nothing to compare with, 
nevertheless our plots show the likely systematic errors resulting 
from the model approximations. 

In Fig. [4] we show (on a logarithmic graph) the integral dis- 
tribution of the ratios of the fitted values of the deflection angle 
parameter to its original values used for the same lenses in simula- 
tions. Similarly Fig.[5]shows the distribution of the ratios between 
fitted and original virial masses. 



4 DISCUSSION AND CONCLUSIONS 

In simulations we replace each halo taken from the Millennium 
Simulation by a specific kind of matter distribution with the same 
mass and size, acting as non-singular isothermal ellipsoid with fi- 
nite mass. We believe that this particular choice of lens model has 
little impact on our results, because our aim is to investigate the 
influence of the environment on the lens action, not the action it- 
self. In particular a different choice of the relative size of the lens 
core (we always use r\jri = 0.01) may change the properties of 
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Figure 3. The integral probability distribution of the external shear value y 
for acceptable models. The results for models neglecting flux ratios (<x m 
1 1' ) ; in Case 1 are shown as solid lines (blue), and in Case 2 - as dashed 
lines (red). For comparison the results for models taking into account flux 
ratios (o~ m = 0.1) in Case 1 are shown as dotted lines (black). The thick 
lines correspond to lower accuracy (o~i = 0.01, Ol = 0.1) and the thin lines 
- to models of higher accuracy (o~i = 0.0025, o~l = 0.025). 
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Figure 4. The integral probability distribution showing the errors in the fit- 
ted values of characteristics deflection angle cio. Conventions follow Fig. [3] 
do is the deflection angle of the lens used in simulations, djj 1 is its fitted 
value. 



Figure 5. The integral probability distribution showing the errors in the 
fitted values of the lens virial mass. Line/colour conventions follow Fig. [5] 
M v j r is the virial mass of the lens used in simulations. M, h '. is its fitted value. 



the radial caustic of each lens, but has no influence on the action of 
matter outside it. 

We consider lenses with non-singular matter distribution so 
the number of images should be odd in all cases. In our calculations 
(which start from an approximate solution of lens equation) the so- 
lutions with four images are rare, 40 times less frequent than the 
five-image cases. The brightest images in our four-image configura- 
tions are on average brighter by ~ 1 .5 mag as compared to the bright- 
est images in five-image cases. Also the flux ratio of the brightest 
image to the next one is substantially higher (by ~ 0.6 mag ) in the 
four-image configurations. This suggests that our four-image so- 
lutions correspond to source positions close to the caustics, which 
produce close image pairs unresolvable by our approximate method 
of image finding. Since four-image configurations are rare, they 
have little impact on the statistics of the acceptability of fits. 

The matter belonging to gravitationally bound haloes makes 
only a part of the whole matter distribution. The distribution of the 
unbound component we are using is known with limited resolution 
of ^ 2 /j-'Mpc in 3D. Since the influence of low density regions 
averages along the lines of sight on distances of several thousands 
megaparsecs, the averaging inside cells of few megaparsecs across 
probably has little meaning. 

Our simulations cover eight small (3' x 3'), randomly chosen 
regions in the sky. In each of the regions we choose ten lenses, with 
largest Einstein ring for a source at z ~ 2 according to approximate, 
SIS lens formula. These are the strongest lenses in the region, giving 
the largest image separations and having the largest area closed by 
caustics in the source plane. That does not mean that every lensed 
source must be related to one of the lenses investigated, but for its 
random position in the region of interest this would happen in a vast 
majority of cases. 

Our numerical experiment shows that matter close to the line 
of sight and matter in the close vicinity of the lens do have an 
important influence on its lensing properties. Our statement is of 
statistical nature: we have shown than in large part of image con- 
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figurations obtained as a result of simulations taking into account 
the presence of matter concentrations near the lens it is impossible 
to obtain acceptable models neglecting their influence. The prob- 
lem is more severe in case when simulations take into account not 
only lens close neighbours but also matter inhomogeneities along 
the beam. It shows that both close neighbours and matter along 
the line o f sight are important. Our resu lt s are in general agree- 
ment w ith lWambsganss et al.|j2004ll2005b : lb'Aloisio & NataraiarJ 
J201 ]]), but we address different aspect of the problem. 

The quantitative characteristics of the effect depends on the 
required accuracy of the simplified modeling. It certainly depends 
on the source redshift as well, but we have not pursued this aspect 
of the problem, limiting ourselves to sources at z s ~ 2. 

Some examples of the expected rate of success, when us- 
ing simplified lens model to fit realistic image configurations, are 
shown in Tables. 1 and 2. We assume that every position of the 
source has the same probability. The selection of a sample of mul- 
tiple image systems may also depend on the lens amplification, im- 
age configuration and other secondary issues, which are difficult to 
simulate and are neglected. 

Our results in Tables. 1 and 2 take into account the fifth image, 
if present. In real cases the fifth images are not observed or very 
faint and difficult to use in modeling. We have performed some ex- 
tra calculations neglecting the fifth images altogether. Of course it 
is always easier to model a less complicated configuration of im- 
ages, so the rate of success is higher in this case. For instance the 
first column of Table. 1 would read: 0.16, 0.48, 0.9996 (instead of 
0.08, 0.38, 0.9998) and the third column: 0.64, 0.68, 0.9995 (in- 
stead of 0.43, 0.59, 0.9997). The changes are substantial, but do 
not qualitatively affect our conclusion. 

Our calculations give typical values of the external shear, 
which is used by the acceptable simplified models. When mod- 
eling configurations resulting from the full approach (taking into 
account all matter inhomogeneities along the rays) we obtain me- 
dian shear value y ~ 0.09 (compare Fig. [3}. For the single plane 
approach (which neglects lens foreground and background objects) 
we get y~ 0.03. Both numbers weakly depend on the required ac- 
curacy of modeling and do not change when the fifth image is ne- 
glected. Comparison of these two numbers and the inspection of 
Table. 1 suggest that shear produced by the matter along the beam 
and the matter in the lens vicinity cannot be neglected in majority 
of cases. It also shows that contribution to the shear by background 
and foreground objects is significant. 

We have also compared the fitted values of lens parameters 
with the parameters of the lens used in the simulation. If the flux 
ratios are not modeled, the fitted parameters may substantially dif- 
fer from their original values. For instance the characteristic bend- 
ing angle of the lens oco may differ by a factor greater than ~ 1.2 
in 50% of cases (compare Fig. [4j . This would imply similar errors 
in estimates of mass distribution in the central part of the lens. For 
the virial mass (M v ; r oco>™-) the discrepancies may exceed factor 
~ 1.5 in 50% of cases (compare Fig. [5). The models reproducing 
also the flux ratios give much better accuracy of the intrinsic lens 
parameters fitting (dotted lines in Fig. [4] [5}. 
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